Phase equilibria and fractionation in a polydisperse fluid 
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We describe how Monte Carlo simulation within the grand canonical ensemble can be applied to 
the study of phase behaviour in polydisperse fluids. Attention is focused on the case of fixed poly- 
dispersity in which the form of the 'parent' density distribution p°(cr) of the polydisperse attribute 
a is prescribed. Recently proposed computational methods facilitate determination of the chemical 
potential distribution conjugate to p°{o). By additionally incorporating extended sampling tech- 
niques within this approach, the compositions of coexisting ('daughter') phases can be obtained and 
fractionation effects quantified. As a case study, we investigate the liquid-vapor phase equilibria of 
a size-disperse Lennard- Jones fluid exhibiting a large (8 = 40%) degree of polydispersity. Cloud and 
shadow curves are obtained, the latter of which exhibit a high degree of fractionation with respect 
to the parent. Additionally, we observe considerable broadening of the coexistence region relative 
to the monodisperse limit. 

I. INTRODUCTION 

Many complex fluids, whether natural or synthetic in origin, comprise mixtures of similar rather than identical 
constituents. Examples are to be found in a host of soft matter materials. For instance, a colloidal dispersion may 
contain particles which exhibit a range of sizes, surface charge or chemical character; while many synthetic materials 
contain macromolecules having a range of chain lengths. This dependence of particle properties on one or more 
continuous parameters is termed polydispersity. It affects the performance of materials in applications ranging from 
foodstuffs to polymer processing [1]. 

In describing polydisperse systems, it is usual to label the polydisperse attribute by a continuous variable a. The 
state of the system is then specified by a density distribution p(a) measuring the number density of particles of 
each a. Certain systems such as micelles, may exhibit variable polydispersity in which the form of p{o~) depends on 
the prevailing chemical and thermodynamic conditions. Others such as colloids and polymers exhibit so-called fixed 
polydispersity because p(a) is set by the synthesis of the fluid. In this letter we shall focus on the latter case. 

Polydisperse fluids differ from their monodisperse counterparts in a variety of aspects. Principal among these 
is the much richer character of their phase behaviour [2]. This richness is traceable to fractionation effects. At 
phase coexistence, particles of each a may partition themselves unevenly between two (or more) coexisting 'daughter' 
phases as long as-due to particle conservation-the overall composition p°(a) of the 'parent' phase is maintained. This 
partitioning alters the character of phase diagrams. For example, the conventional liquid-gas binodal of a monodisperse 
system (which connects the ends of tie-lines in a density-temperature diagram) splits into a 'cloud' and a 'shadow' 
curve. These give, respectively, the density at which phase coexistence first occurs and the density of the incipient 
phase; the curves do not coincide because the shadow phase in general differs in composition from the parent. Only 
recently has experimental work started to elucidate in a systematic fashion the generic consequences of fractionation 
for phase coexistence properties [3-5]. 

Computational solutions for dealing with polydispersity have generally focused on the semi-grand canonical ensemble 
(SGCE). Within this framework, the instantaneous form of p(a) is permitted to fluctuate under the control of a 
distribution of chemical potential differences pd(&), subject to a fixed overall number of particles. Use of such an 
approach is attractive because it permits the sampling of many different realizations of the ensemble of particle sizes, 
thereby ameliorating finite-size effects. For the investigation of phase coexistence, the SGCE has been combined with 
Gibbs-Duhem integration [6, 7] and Gibbs ensemble simulations [8, 9]. However, these studies were restricted to the 
case of variable polydispersity; no attempts were made to target a specific form of p{o~) or determine cloud and shadow 
curves. 

The computational difficulties associated with tackling fixed polydispersity are potentially quite severe: one needs 
to determine that form of the chemical potential distribution p(o~) for which the ensemble averaged composition dis- 
tribution p(o~) matches the target i.e. the prescribed parent p°(o~). Unfortunately, the chemical potential distribution 
is unavailable a priori, it being an unknown functional (i.e. p(cr) = p[p°(a)}), of the parent. Recently however, tech- 
niques have been developed that efficiently overcome this difficulty within the framework of a full grand canonical 
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ensemble (GCE) [10-12]. The latter is particular well suited to the study of fluid phase transitions due to the fluctu- 
ating overall particle number. In this letter we demonstrate that when combined with extended sampling methods, 
the new techniques facilitate the detailed and efficient study of phase behaviour in fluids of fixed polydispersity. 

II. COMPUTATIONAL ASPECTS 

The grand canonical ensemble Monte Carlo algorithm we employ has been described in ref. [10] and invokes four 
types of operation: particle displacements, deletions, insertions, and resizing. The polydisperse attribute a is itself 
represented as a strictly continuous variable, subject to some upper bound <r c . However, observables such as the 
instantaneous composition distribution p(a) are accumulated in the form of a histogram by discretising the a domain 
into a prescribed number of bins. This discretisation also applies to the chemical potential distribution /u(«r), i.e. all 
particles whose a values is encompassed by the same bin are subject to an identical chemical potential. 

The form of the ensemble averaged composition distribution p(a) is controlled by /x(cr), via its role in the acceptance 
probabilities for particle transfers and resizing moves. For fixed polydispersity one wishes to match p(a) to the desired 
parent distribution. The latter can be written as 

p°(a) = n /(a), (1) 

where no = N/V is the overall particle number density, while /(c) is a prescribed normalized shape function. Since 
p°(cr) may vary only in terms of its scale no, the system is constrained to traverse a dilution line in the full phase 
space of possible compositions. The task is then to determine, as a function of temperature, the form of p(a) along 
the dilution line. More specifically, we seek the intersection of the dilution line with a coexistence region. Recently 
developed simulation techniques facilitate this, as we now summarize. 

The non-equilibrium potential refinement (NEPR) scheme [12] permits the efficient iterative determination of 
p[p°(a),T], from a single simulation run, and without the need for an initial guess of its form. To achieve this, 
the method continually updates p(o~) in such as way as to minimize the deviation of the instantaneous density dis- 
tribution p(a) from the target form (i.e. the parent). However, tuning p(a) in this manner clearly violates detailed 
balance. To counter this, successive iterations reduce the degree of modification applied to u(er), thereby driving the 
system towards equilibrium and ultimately yielding the equilibrium form of p[p°(o~), T]. 

For the purpose of exploring phase diagrams, Histogram Extrapolation (HE) techniques have proved invaluable 
[13]. In the present context, their use permits histogram of observables accumulated at one p(a) to be reweighted to 
estimate observables at some other p{o~). In ref. [10] we have shown how HE can be combined with a minimization 
scheme, to track a dilution line in a stepwise fashion. We shall deploy this approach again in the present study. 

Simulation studies of phase coexistence present distinctive challenges. Principal among these is the large free 
energy (surface tension) barrier separating the coexisting phases. In order to accurately locate coexistence points, a 
sampling scheme must be utilized which enables this barrier to be surmounted [14]. One such scheme is multicanonical 
preweighting, which utilizes a weight function in the MC acceptance probabilities, in order to encourage the simulation 
to sample the interfacial configurations of low probability [15]. At a given coexistence state point, the requisite weight 
function takes the form of an approximation to the inverse of the distribution of the fluctuating number density, p(n) , 
with n = N/V. While specialized techniques allow determination of p(n) from scratch, in situations where one wishes 
to track a fluid-fluid phase boundary, prior determination of a weight function is unnecessary provided one commences 
from the vicinity of the critical point where the barrier to inter-phase crossings is small. Data accumulated here can 
be used (together with HE) to provide estimates of suitable multicanonical weight functions at lower temperatures 
[16] where the barrier height is greater. 

III. MODEL 

The techniques outlined above have been deployed to obtain the liquid-vapor coexistence properties of a fluid of 
particles interacting via a pairwise potential of the Lennard- Jones (LJ) form: 



(2) 

with cry = (<7j + o-j)/2. A cutoff was applied to this potential for particle separations > 2.5(7^. 



u ( r i 3 ,o-ij) = 4e 
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The polydispersity enters solely through the distribution of diameters crj. Our algorithm finds that n(cr) for which 
the diameters are distributed according to /(<r) given a choice for n (cf. cq. 1). We have assigned f(a) the Schulz 
form: 



z+l 



a exp 



(3) 



Here a = 1 is the average particle diameter, while z is a width parameter, the value of which was set to z = 5, 
corresponding to a high degree of polydispersity (standard deviation of /(c)), S = 40%. Additionally, for convenience, 
f(a) was truncated at a c — 3.0. Histograms of observables were formed by discretising the permitted range < a < a c 
into 120 bins. 



IV. SIMULATION STRATEGY AND RESULTS 



The strategy employed for mapping the liquid-vapor coexistence curve of our model was as follows. In order to 
bootstrap the dilution line tracking procedure, the NEPR method [12] was employed to determine \i{a) for a gas phase 
state point on the dilution line at a moderately low temperature. Starting from this point, the dilution line was then 
followed towards increasing density (with the aid of HE) until the gas spontaneously liquefied. Having estimated the 
location of a coexistence state point in this manner, the temperature was increased in steps (whilst remaining on the 
dilution line) until the density difference between the gas and the spontaneously formed liquid vanished, signalling the 
proximity of the critical point. Finite-size scaling methods [16] were then used to home in on the critical parameters. 




FIG. 1: The measured phase diagram of the size-disperse LJ fluid, (a) The n — T representation, (b) The r\ — T representation. 
In both cases, the phase diagram of the monodisperse limit (broken line) is shown for comparison. Critical points (determined 
by finite-size scaling) are shown as crosses. Statistical errors do not exceed the symbol sizes. 



Having located the critical point, a detailed mapping of the cloud and shadow curves was performed for a large 
simulation box of volume V = 11390ct 3 . Attention was focused on the distribution of the fluctuating overall number 
density, p(n). The gas phase cloud point (incipient liquid phase) corresponds to the situation where p(n) is bimodal, 
but with vanishingly small weight in the liquid peak. Under these conditions, the position of the low density gas peak 
provides an estimate of the gas phase cloud density, while that of the liquid peak gives the gas phase shadow density. 
The converse is true for the liquid phase cloud point and its shadow. Determining the cloud and shadow points as a 
function of temperature yields the cloud and shadow curves. We have tracked the gas and liquid cloud curves (and 
their shadows) in a stepwise fashion downwards in temperature from the critical point. Histogram extrapolation was 
employed to negotiate each temperature step, yielding estimates for both the form of /i(a) on the cloud curve at the 
next temperature, and the requisite multicanonical weight function. 

The resulting phase diagram is shown in fig. 1(a), together with that of the monodisperse LJ fluid, determined in an 
earlier study [16]. It should be pointed out that while the positions of the peaks in p(n) provide an accurate estimate 
of cloud and shadow points at low temperatures, this breaks down near the critical point due to finite-size effects [16]. 
Thus a naive extrapolation of our curves to their intersection point will tend to overestimate the critical temperature. 
However, our independent determination of the critical point using finite-size scaling methods (as indicated in fig. 1(a)) 
is considerably more accurate. 
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The results of fig. 1 (a) show a stark separation of the cloud and shadow curves in the n — T plane. Furthermore, the 
whole phase diagram is considerably shifted with respect to that of the monodispcrsc fluid. Specifically, one observes 
that the critical point occurs at a considerably higher temperature than in the monodisperse limit. This particular 
finding contrasts with that of a previous theoretical study of a size-disperse van-der Waals fluid [17], which predicts 
a suppression of the critical temperature with respect to the monodisperse limit. 
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FIG. 2: The measured form of the chemical potential distribution n(a) at the gas and liquid phase cloud points for T — 0.91T C . 
Statistical errors do not exceed the symbol sizes. 

The standard order of cloud and shadow curves with increasing density is cloud-shadow-cloud-shadow. By contrast, 
the pattern apparent in fig. 1(a) is cloud-shadow-shadow-cloud. Interestingly, however, the order reverts to the 
standard pattern if one plots the data in terms of the volume fraction r\ = (n/6) J daa 3 p(a), rather than the number 
density, as shown in fig. 1(b). Moreover, one sees that in the r\ — T representation the differences between cloud 
and shadow phase properties become much less pronounced. In particular, while the critical number density of the 
polydisperse fluid is considerably less than its value in the monodisperse limit, the critical volume fraction for the 
mono- and polydisperse fluid agree to within error. However, irrespective of the choice of data representation, we 
observe that for our model the critical point occurs very close to the top of the coexistence curve. No clear evidence 
was discernible for distinct cloud and shadow points, at or above T c . 

Notwithstanding these intriguing findings, not all differences between gas and liquid cloud points at a given temper- 
ature can be camouflaged by a simple change of variable. At temperatures significantly below criticality, we observe 
dramatic broadening of the coexistence curve in the space of /x(cr). This is shown in fig. 2 which presents the form of 
fi(a) at the respective cloud points for the lowest temperature studied, T = 0.91T C . Such broadening does not occur 
in monodisperse systems (coexistence occurs at a single value of the chemical potential, not a range of values). The 
effect is surprisingly large, even given the high degree of polydispersity of the parent (8 — 40%). Indeed, in simulation 
terms, the respective cloud points are so far separated in phase space that to connect them directly (via a route 
crossing the phase boundary) required a dozen overlapping simulations-twice as many as were required to connect 
the cloud point to the critical point at this temperature. We remark in passing that similar aspects of coexistence 
curve broadening have recently been analyzed within the context of Landau theory by Rascon and Cates [18]. 

Finally in this section we present the normalized daughter phase distributions at the gas and liquid cloud points for 
T = 0.91T C . The data show that at the gas phase cloud point, larger particles preferentially occupy the shadow phase. 
Conversely at the liquid phase cloud point, there is a predominance of smaller particles in the shadow phase. Clearly 
the scale of these fractionation effects is considerable: the polydispersity of the gas phase shadow at this temperature 
is close to 50%, while that of the liquid phase shadow is « 33%, to be compared with a parent polydispersity of 40%. 

V. DISCUSSION AND CONCLUSIONS 

In summary, we have demonstrated how extended sampling grand canonical simulations can be combined with 
histogram extrapolation methods and a new non-equilibrium potential refinement scheme to accurately determine the 
phase behaviour of a polydisperse fluid. The results show that in contrast to existing theoretical predictions, the critical 
temperature of the polydisperse system exceeds that of its monodisperse counterpart. As regards the sub-critical 
region, we find that the relative order of cloud and shadow curves changes depending on whether the data is represented 
in terms of the overall number density or the volume fraction. Additionally, we observe considerable polydispersity- 
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FIG. 3: The normalized form of the particle size distribution at the gas and liquid phase shadow points at T = 0.91T C . The 
average particle diameter in the gas phase shadow is a — 1.167(3), while that in the liquid-phase shadow is a — 0.863(2). Also 
shown for comparison is the parent shape function f(a), corresponding to the respective cloud phase distributions. 

induced broadening of the coexistence region: at a given temperature, the cloud points of the respective phases (which 
coincide in a monodisperse system) are widely separated in terms of their chemical potential distributions. The scale 
of this effect is mirrored in the disparate forms of the shadow phase daughter distributions. 

In a future publication [19], we will present further simulation results for the magnitude of critical point shifts 
as a function of the width of the governing size distribution f(cr). These results will be compared with those of a 
moment based theoretical analysis [20, 21] of an improved model free energy for the size disperse van der Waals fluid. 
The latter correctly captures the sign of polydispersity-induced critical point shifts and provides insights into the 
deficiencies of previous approaches. 
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